Here I will explore how DEPOSITED AND SUSPENDED sediment affects the PHOTOSYNTHESIS of tropical, scleractinean coral ADULTS. The database to be used is comprised of data extracted from studies deemed relevant during the systematic literature review process.

The Dataset

Deposited Sediment

Specifically, we will explore the results of 438 datapoints from 15 studies conducted in 3 oceans on 36 species within 23 genera.

##     n
## 1 416

##     Ref   n
## 1  DS34 157
## 2  DS04  58
## 3  DS68  48
## 4  DS16  46
## 5  DS17  24
## 6  DS40  20
## 7  DS10  16
## 8  DS03  12
## 9  DS43  12
## 10 DS20   8
## 11 DS12   7
## 12 DS24   4
## 13 DS28   4

##                            Ref_name   n
## 1              Weber et al. (2006)  157
## 2           Duckworth et al. (2017)  58
## 3              Flores et al. (2012)  48
## 4      Philipp and Fabricius (2003)  46
## 5                     Piniak (2007)  24
## 6  Abdel-Salam (1989) PhD Chapter 3  20
## 7              Junjie et al. (2014)  16
## 8    Bessell-Browne et al. (2017a)   12
## 9         Peters and Pilson (1985)   12
## 10          Riegl and Branch (1995)   8
## 11             Loiola et al. (2013)   7
## 12           Sheridan et al. (2014)   4
## 13       Sofonia and Anthony (2008)   4

##          Ocean   n
## 1      Pacific 349
## 2     Atlantic  39
## 3 Indo-Pacific  16
## 4       Indian  12

## Selecting by n

##                           Gsp   n
## 1          Acropora_millepora  28
## 2          Astrangia_poculata  12
## 3        Galaxea_fascicularis  11
## 4       Goniopora_somaliensis   8
## 5  Montipora_aequituberculata  24
## 6          Montipora_capitata  12
## 7       Montipora_capricornis  20
## 8       Montipora_peltiformis 167
## 9        Porites_lobata/lutea  39
## 10      Turbinaria_reniformis  25

## Selecting by n

##    Updated_Genus   n
## 1      Montipora 236
## 2        Porites  43
## 3     Turbinaria  32
## 4       Acropora  30
## 5      Astrangia  12
## 6        Galaxea  11
## 7      Goniopora   8
## 8    Mussismilia   7
## 9      Orbicella   6
## 10     Ctenactis   3
## 11    Echinopora   3
## 12      Merulina   3
## 13    Pachyseris   3
## 14      Pectinia   3

Suspended Sediment

Specifically, we will explore the results of 251 datapoints from 7 studies conducted in 3 oceans on 9 species within 6 genera (shown below only the top 10 species and genera, by number of datapoints).

##     n
## 1 251

##     Ref   n
## 1  SS06 146
## 2  SS08  46
## 3 SS17b  20
## 4 SS17a  17
## 5  SS05  12
## 6  SS04   9
## 7  SS25   1

##                     Ref_name   n
## 1         Browne et al. 2015 146
## 2         Flores et al. 2012  46
## 3            Liu et al. 2015  37
## 4         Browne et al. 2014  12
## 5 Bessell-Browne et al. 2017   9
## 6          Te 2001 Chapter 6   1

##                           Ocean   n
## 1                        Indian 146
## 2                       Pacific  93
## 3 Indo-Pacific (e.g. Singapore)  12

##                          Gsp  n
## 1         Platygyra_sinensis 54
## 2          Merulina_ampliata 53
## 3        Pachyseris_speciosa 51
## 4          Acropora_muricata 37
## 5 Montipora_aequituberculata 30
## 6         Acropora_millepora 19
## 7      Montipora_capricornis  3
## 8      Porites_lutea, lobata  3
## 9        Montipora_verrucosa  1

##   Updated_Genus  n
## 1      Acropora 56
## 2     Platygyra 54
## 3      Merulina 53
## 4    Pachyseris 51
## 5     Montipora 34
## 6       Porites  3

Exploratory plots

First let’s explore all data from all species for which there is binary data about the presence of a reduction in ‘photosynthetic efficiency’ or ‘P/R ratio’ as a result of exposure to deposited sediment.

DEFINITIONS:

‘Photosynthetic efficiency’ is defined as either maximum or effective quantum yield, reported as Fv/Fm or deltaF/F’m, respectively.
‘P/R ratio’ is defined as the ratio between photosynthetic and respiration rates.

Threshold Analyses with Binary Data

Now let’s calculate the thresholds for each category of photosynthesis, based on the binary data explored above.

DEFINITIONS:

‘LOAEL’, or the ‘lowest observed adverse effect level’ is defined as the lowest dose at which there was an observed adverse effect.
‘NOAEL’, or the ‘no observed adverse effect level’ is defined as the highest dose at which there was NOT an observed adverse effect.
‘Adverse effect’ is defined here as any response of a coral individual, colony, or experimental treatment group that may negatively affect the coral’s fitness and/or survival. These adverse effects may include sublethal physiological changes (e.g., significantly reduced growth or photosynthetic rates when compared to coral in ambient/control conditions), bleaching/tissue loss, and mortality. This definition of adverse effect is independent of its magnitude, so while the affect may have potentially reduced fitness, its magnitude may be sufficiently small that the fitness effect is not measureable.

DEPOSITED SEDIMENT

Photosynthetic efficiency

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1    25

##   NOAEL
## 1    25

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1   0.5

##   NOAEL
## 1   0.5

P/R ratio

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1  26.4

##   NOAEL
## 1  26.4

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1     2

##   NOAEL
## 1     2

SUSPENDED SEDIMENT

Photosynthetic efficiency

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1  35.8

##   NOAEL
## 1  35.8

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1    56

##   NOAEL
## 1    56

P/R Ratio

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1  35.8

##   NOAEL
## 1  35.8

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1 0.083

##   NOAEL
## 1 0.083

Logistic Regression Model Fitting

DEPOSITED SEDIMENT

Photosynthetic efficiency

##     n
## 1 249

## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular

## Data: photo_effDS2
## Models:
## modDS6: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS6:     (1 | Gsp)
## modDS9: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS9:     (1 | Ref)
## modDS12: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS12:     (1 | Gsp) + (1 | Ref)
##         npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## modDS6     5 327.12 344.71 -158.56   317.12                          
## modDS9     5 314.54 332.13 -152.27   304.54 12.5790  0     <2e-16 ***
## modDS12    6 314.16 335.26 -151.08   302.16  2.3866  1     0.1224    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Data: photo_effDS2
## Models:
## modDS7: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Ref)
## modDS8: Binary_reduced_FvFm ~ Sed_level_stand_mg + Sed_exposure_days + 
## modDS8:     (1 | Ref)
## modDS9: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS9:     (1 | Ref)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## modDS7    3 312.94 323.49 -153.47   306.94                     
## modDS8    4 313.04 327.11 -152.52   305.04 1.8973  1     0.1684
## modDS9    5 314.54 332.13 -152.27   304.54 0.4953  1     0.4816

## Data: photo_effDS2
## Models:
## modDS10: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## modDS11: Binary_reduced_FvFm ~ Sed_level_stand_mg + Sed_exposure_days + 
## modDS11:     (1 | Gsp) + (1 | Ref)
## modDS12: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS12:     (1 | Gsp) + (1 | Ref)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## modDS10    4 313.93 328.00 -152.97   305.93                       
## modDS11    5 312.74 330.33 -151.37   302.74 3.1894  1    0.07412 .
## modDS12    6 314.16 335.26 -151.08   302.16 0.5867  1    0.44371  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Data: photo_effDS2
## Models:
## modDS4: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp)
## modDS7: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Ref)
## modDS10: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## modDS11: Binary_reduced_FvFm ~ Sed_level_stand_mg + Sed_exposure_days + 
## modDS11:     (1 | Gsp) + (1 | Ref)
##         npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## modDS4     3 328.68 339.23 -161.34   322.68                          
## modDS7     3 312.94 323.49 -153.47   306.94 15.7408  0    < 2e-16 ***
## modDS10    4 313.93 328.00 -152.97   305.93  1.0031  1    0.31656    
## modDS11    5 312.74 330.33 -151.37   302.74  3.1894  1    0.07412 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Ref)
##    Data: photo_effDS2
## 
##      AIC      BIC   logLik deviance df.resid 
##    312.9    323.5   -153.5    306.9      246 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8252 -0.8127 -0.3036  1.0864  3.3269 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Ref    (Intercept) 2.737    1.655   
## Number of obs: 249, groups:  Ref, 9
## 
## Fixed effects:
##                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)        -1.0456482  0.6970200  -1.500    0.134
## Sed_level_stand_mg -0.0008905  0.0022947  -0.388    0.698

Effect estimates

## [1] 0.6706827

##    
## p     0   1
##   0 148  76
##   1   6  19

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Area under the curve: 0.7535

##                     R2m       R2c
## theoretical 0.001289846 0.4548780
## delta       0.001055897 0.3723734

## $Ref

##  Groups Name        Std.Dev.
##  Ref    (Intercept) 1.6545

## [1] 5.230593

There is no significant relationship between deposited sediment exposure and the odds of reduced photosynthetic efficiency (GLMM z = -0.388, p = 0.698).

Prediction figure

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDS10 <- predict(modDS10, newdata=photo_effDS2, type="response")
photo_effDS3 <- cbind(photo_effDS2, pred_modDS10)

photo_effDS_plot <- ggplot(data = photo_effDS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_reduced_FvFm)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Reduced photosynthetic efficiency\ndue to sediment exposure?") +
    scale_x_log10(limits=c(0.01,1000),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modDS10), inherit.aes=FALSE) +
    theme_bw()

photo_effDS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

jvalues <- with(photo_effDS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    photo_effDS2$Sed_level_stand_mg <- j
    predict(modDS10, newdata = photo_effDS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
photo_effDS_plot2 <- ggplot() +
    geom_point(data = photo_effDS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_reduced_FvFm,
      color=Ref)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Predicted probability of reduced photosynthetic\nefficiency due to sediment exposure",
         color = "Binary data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(0.01,1000), breaks=c(0.01,0.1,1,10,100,1000,loaelDS1),
                  label=c("0.01","0.1","1","10","100","1000",round(loaelDS1,digits=1))) +
    geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelDS1, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

photo_effDS_plot2

P/R Ratio

##    n
## 1 25

##    Ref n
## 1 DS10 8
## 2 DS24 2
## 3 DS40 9
## 4 DS43 6

##     Ref                   Gsp
## 1  DS10  Galaxea_fascicularis
## 5  DS10 Goniopora_somaliensis
## 9  DS24      Montipora_patula
## 11 DS40   Agaricia_agaricites
## 12 DS40       Porites_porites
## 13 DS40    Porites_astreoides
## 14 DS40     Diploria_strigosa
## 15 DS40     Isophyllia_rigida
## 16 DS40   Orbicella_annularis
## 17 DS40  Acropora_cervicornis
## 18 DS40  Meandrina_meandrites
## 20 DS43    Astrangia_poculata

## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular

## Data: PR_ratioDS2
## Models:
## modDSb6: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb6:     (1 | Gsp)
## modDSb9: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb9:     (1 | Ref)
## modDSb12: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb12:     (1 | Gsp) + (1 | Ref)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modDSb6     5 16.530 22.625 -3.2652   6.5304                         
## modDSb9     5 13.379 19.474 -1.6895   3.3791 3.1513  0     <2e-16 ***
## modDSb12    6 15.379 22.693 -1.6897   3.3793 0.0000  1          1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Data: PR_ratioDS2
## Models:
## modDSb7: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Ref)
## modDSb8: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + Sed_exposure_days + 
## modDSb8:     (1 | Ref)
## modDSb9: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb9:     (1 | Ref)
##         npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## modDSb7    3 25.440 29.096 -9.7198  19.4396                          
## modDSb8    4 12.511 17.386 -2.2554   4.5108 14.9288  1  0.0001116 ***
## modDSb9    5 13.379 19.474 -1.6895   3.3791  1.1317  1  0.2874148    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Data: PR_ratioDS2
## Models:
## modDSb4: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Gsp)
## modDSb7: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Ref)
## modDSb10: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modDSb10:     Ref)
##          npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)    
## modDSb4     3 32.834 36.491 -13.4171   26.834                         
## modDSb7     3 25.440 29.096  -9.7198   19.440 7.3946  0     <2e-16 ***
## modDSb10    4 27.440 32.315  -9.7198   19.440 0.0000  1     0.9994    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + Sed_exposure_days +  
##     (1 | Ref)
##    Data: PR_ratioDS2
## 
##      AIC      BIC   logLik deviance df.resid 
##     12.5     17.4     -2.3      4.5       21 
## 
## Scaled residuals: 
##        Min         1Q     Median         3Q        Max 
## -0.0005026 -0.0000769  0.0000000  0.0111328  0.0111328 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Ref    (Intercept) 78355    279.9   
## Number of obs: 25, groups:  Ref, 4
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -76.1954    20.5559  -3.707 0.000210 ***
## Sed_level_stand_mg  -0.0217     0.1567  -0.138 0.889872    
## Sed_exposure_days    4.3993     1.1938   3.685 0.000228 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 1

##    
## p    0  1
##   0 12  0
##   1  0 13

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Area under the curve: 1

##                    R2m       R2c
## theoretical 0.01728004 0.9999587
## delta       0.01725835 0.9987037

## $Ref

##  Groups Name        Std.Dev.
##  Ref    (Intercept) 279.92

## [1] 3.694971e+121

##                             Est           LL           UL
## (Intercept)        1.155860e-23 8.600357e-36 1.553438e-11
## Sed_level_stand_mg 9.850746e-01 7.962031e-01 1.218749e+00
## Sed_exposure_days  2.110169e+01 4.168496e+00 1.068206e+02

For every doubling in exposure duration of deposited sediment, the odds of reduced P/R ratio increases by 21.1 times (95% CI 4.2, 106.8; GLMM z = 3.685, p = 0.0002), after accounting for exposure concentration. There is no significant relationship between sediment exposure concentration and the odds of reduced P/R ratio (GLMM z = -0.138, p = 0.890). However, all the variability in the model that describes this relationship is attributable to study (R2 = 0.017 with fixed effects only, R2 = 1.000 with fixed effects + random effect of study).

Prediction figure

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDSb8 <- predict(modDSb8, newdata=PR_ratioDS2, type="response")
PR_ratioDS3 <- cbind(PR_ratioDS2, pred_modDSb8)

PR_ratioDS_plot <- ggplot(data = PR_ratioDS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_reduced_PR_ratio,
      color = Ref)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Reduced P/R ratio due to sediment exposure?",
         color = "Study") +
    scale_x_continuous(limits=c(0,200)) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modDSb8), inherit.aes=FALSE) +
    theme_bw()

PR_ratioDS_plot

By exposure concentration

jvalues <- with(PR_ratioDS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    PR_ratioDS2$Sed_level_stand_mg <- j
    predict(modDSb8, newdata = PR_ratioDS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
PR_ratioDS_plot2 <- ggplot() +
    geom_point(data = PR_ratioDS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_reduced_PR_ratio,
      color = Ref)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Predicted probability of reduced\nP/R ratio due to sediment exposure",
         color = "Binary data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_continuous(limits=c(0,200), breaks=c(0,50,100,150,200,loaelDS2),
                       labels=c("0","50","100","150","200",round(loaelDS2,digits=1))) +
    geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelDS2, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

PR_ratioDS_plot2

By exposure duration

jvalues <- with(PR_ratioDS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    PR_ratioDS2$Sed_exposure_days <- j
    predict(modDSb8, newdata = PR_ratioDS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_dur)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
PR_ratioDS_plot2 <- ggplot() +
    geom_point(data = PR_ratioDS2, mapping = aes(
      x = Sed_exposure_days, 
      y = Binary_reduced_PR_ratio,
      color = Ref)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Predicted probability of reduced\nP/R ratio due to sediment exposure",
         color = "Study",
         linetype = "Predicted\nProbability") +
    scale_x_continuous(limits=c(0,30)) +
    geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability, 
                                  linetype = "solid")) +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

PR_ratioDS_plot2

SUSPENDED SEDIMENT

Photosynthetic Efficiency

##     n
## 1 180

## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular

## Data: photo_effSS2
## Models:
## modSS6: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS6:     (1 | Gsp)
## modSS9: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS9:     (1 | Ref)
## modSS12: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS12:     (1 | Gsp) + (1 | Ref)
## modSS15: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS15:     (1 | Ref_name/Ref)
## modSS18: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS18:     (1 | Gsp) + (1 | Ref_name)
## modSS21: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS21:     (1 | Gsp) + (1 | Ref_name/Ref)
##         npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)    
## modSS6     5 29.924 45.889  -9.9622   19.924                         
## modSS9     5 36.257 52.222 -13.1285   26.257 0.0000  0    1.00000    
## modSS12    6 31.924 51.082  -9.9622   19.924 6.3326  1    0.01185 *  
## modSS15    6 38.256 57.414 -13.1281   26.256 0.0000  0    1.00000    
## modSS18    6 31.924 51.082  -9.9622   19.924 6.3316  0    < 2e-16 ***
## modSS21    7 33.924 56.275  -9.9622   19.924 0.0000  1    1.00000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Data: photo_effSS2
## Models:
## modSS4: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp)
## modSS5: Binary_reduced_FvFm ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSS5:     (1 | Gsp)
## modSS6: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS6:     (1 | Gsp)
##        npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)   
## modSS4    3 33.228 42.807 -13.6138   27.228                        
## modSS5    4 28.037 40.808 -10.0183   20.037 7.1911  1   0.007327 **
## modSS6    5 29.924 45.889  -9.9622   19.924 0.1121  1   0.737798   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Data: photo_effSS2
## Models:
## modSS4: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp)
## modSS7: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Ref)
## modSS10: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## modSS13: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## modSS16: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref_name)
## modSS19: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref_name/Ref)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modSS4     3 33.228 42.807 -13.614   27.228                         
## modSS7     3 36.149 45.728 -15.075   30.149 0.0000  0     1.0000    
## modSS10    4 35.228 47.999 -13.614   27.228 2.9215  1     0.0874 .  
## modSS13    4 38.149 50.921 -15.075   30.149 0.0000  0     1.0000    
## modSS16    4 35.228 47.999 -13.614   27.228 2.9215  0     <2e-16 ***
## modSS19    5 37.228 53.192 -13.614   27.228 0.0000  1     1.0000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_reduced_FvFm ~ Sed_level_stand_mg + Sed_exposure_days +  
##     (1 | Gsp)
##    Data: photo_effSS2
## 
##      AIC      BIC   logLik deviance df.resid 
##     28.0     40.8    -10.0     20.0      176 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.20089 -0.00726 -0.00048 -0.00015  2.64336 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Gsp    (Intercept) 94.68    9.731   
## Number of obs: 180, groups:  Gsp, 8
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        -21.30065    9.94516  -2.142   0.0322 *
## Sed_level_stand_mg   0.02395    0.02453   0.976   0.3289  
## Sed_exposure_days    0.13636    0.07874   1.732   0.0833 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 0.9777778

##    
## p     0   1
##   0 176   3
##   1   1   0

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Area under the curve: 0.9859

##                    R2m       R2c
## theoretical 0.09587058 0.9696399
## delta       0.08433380 0.8529564

## $Gsp

##  Groups Name        Std.Dev.
##  Gsp    (Intercept) 9.7305

## [1] 16823.41

##                             Est           LL        UL
## (Intercept)        3.871372e-07 5.248454e-13 0.2855607
## Sed_level_stand_mg 1.016741e+00 9.834110e-01 1.0512007
## Sed_exposure_days  1.099130e+00 9.876226e-01 1.2232273

There is no significant relationship between exposure concentration of suspended sediment and the odds of reduced photosynthetic efficiency (GLMM z = 0.976, p = 0.329), though there is suggestive, but non-significant evidence that for every doubling of exposure duration that the odds of reduced photosynthetic efficiency increase by 1.1 times (95% CI 1.0, 1.2; GLMM z = 1.732, p = 0.083).

Prediction figures

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSS5 <- predict(modSS5, newdata=photo_effSS2, type="response")
photo_effSS3 <- cbind(photo_effSS2, pred_modSS5)

photo_effSS_plot <- ggplot(data = photo_effSS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_reduced_FvFm,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Reduced photosynthetic efficiency due to sediment exposure",
         color = "Study") +
    scale_x_log10(limits=c(1,250),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modSS5), inherit.aes=FALSE) +
    theme_bw()

photo_effSS_plot

By sediment exposure concentration

jvalues <- with(photo_effSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    photo_effSS2$Sed_level_stand_mg <- j
    predict(modSS5, newdata = photo_effSS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
photo_effSS_plot2 <- ggplot() +
    geom_point(data = photo_effSS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_reduced_FvFm,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Predicted probability of reduced photosynthetic\nefficiency due to sediment exposure",
         color = "Binary Data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(1,max(photo_effSS2$Sed_level_stand_mg)), 
                  breaks=c(0.01,0.1,1,10,100,1000,loaelSS1),
                  label=c("0.01","0.1","1","10","100","1000",round(loaelSS1,digits = 1))) +
    geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelSS1, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

photo_effSS_plot2

By sediment exposure duration

jvalues <- with(photo_effSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    photo_effSS2$Sed_exposure_days <- j
    predict(modSS5, newdata = photo_effSS2, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_dur)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
photo_effSS_plot3 <- ggplot() +
    geom_point(data = photo_effSS2, mapping = aes(
      x = Sed_exposure_days, 
      y = Binary_reduced_FvFm,
      color = Ref)) +
    labs(x = "Sediment exposure duration (days)", 
         y = "Predicted probability of reduced photosynthetic\nefficiency due to sediment exposure",
         color = "Binary Data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_continuous(limits=c(0,max(photo_effSS2$Sed_exposure_days))) +
    geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability, 
                                  linetype = "solid")) +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

photo_effSS_plot3

P/R Ratio

##    n
## 1 34

## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular

## Data: PR_ratioSS2
## Models:
## modSSb6: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSSb6:     (1 | Gsp)
## modSSb9: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSSb9:     (1 | Ref)
## modSSb12: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSSb12:     (1 | Gsp) + (1 | Ref)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## modSSb6     5 36.328 43.959 -13.164   26.328                     
## modSSb9     5 36.378 44.010 -13.189   26.378 0.0000  0      1.000
## modSSb12    6 38.328 47.486 -13.164   26.328 0.0506  1      0.822

## Data: PR_ratioSS2
## Models:
## modSSb7: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Ref)
## modSSb8: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSSb8:     (1 | Ref)
## modSSb9: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSSb9:     (1 | Ref)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## modSSb7    3 33.577 38.156 -13.788   27.577                     
## modSSb8    4 35.425 41.531 -13.713   27.425 0.1514  1     0.6972
## modSSb9    5 36.378 44.010 -13.189   26.378 1.0468  1     0.3062

## Data: PR_ratioSS2
## Models:
## modSSb4: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Gsp)
## modSSb7: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Ref)
## modSSb10: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modSSb10:     Ref)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## modSSb4     3 33.549 38.128 -13.774   27.549                     
## modSSb7     3 33.577 38.156 -13.788   27.577 0.0000  0     1.0000
## modSSb10    4 35.549 41.654 -13.774   27.549 0.0276  1     0.8682

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Ref)
##    Data: PR_ratioSS2
## 
##      AIC      BIC   logLik deviance df.resid 
##     33.6     38.2    -13.8     27.6       31 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.6581 -0.4039 -0.4039 -0.3360  2.9765 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev. 
##  Ref    (Intercept) 1.256e-16 1.121e-08
## Number of obs: 34, groups:  Ref, 3
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)   
## (Intercept)        -2.414418   0.914690  -2.640   0.0083 **
## Sed_level_stand_mg  0.006506   0.007018   0.927   0.3539   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Gsp)
##    Data: PR_ratioSS2
## 
##      AIC      BIC   logLik deviance df.resid 
##     33.5     38.1    -13.8     27.5       31 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.6967 -0.4137 -0.3646 -0.3060  2.8420 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Gsp    (Intercept) 0.1163   0.341   
## Number of obs: 34, groups:  Gsp, 4
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        -2.468090   1.002623  -2.462   0.0138 *
## Sed_level_stand_mg  0.006610   0.007112   0.929   0.3527  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 0.8529412

##    
## p    0  1
##   0 29  5

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Area under the curve: 0.6655

##                    R2m        R2c
## theoretical 0.04833292 0.08081445
## delta       0.02088213 0.03491570

## $Gsp

##  Groups Name        Std.Dev.
##  Gsp    (Intercept) 0.34096

## [1] 1.4063

There is no significant relationship between exposure concentration of suspended sediment and the odds of a reduced P/R ratio (GLMM z = 0.929, p = 0.353), but the best-fit model to explain this relationship explains only 8% of the variance, indicating that more data is needed.

Prediction figures

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSSb4 <- predict(modSSb4, newdata=PR_ratioSS2, type="response")
PR_ratioSS3 <- cbind(PR_ratioSS2, pred_modSSb4)

PR_ratioSS_plot <- ggplot(data = PR_ratioSS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_reduced_PR_ratio,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Reduced P/R ratio due to sediment exposure",
         color = "Study") +
    scale_x_log10(limits=c(10,1000),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modSSb4), inherit.aes=FALSE) +
    theme_bw()

PR_ratioSS_plot

By sediment exposure concentration

jvalues <- with(PR_ratioSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    PR_ratioSS2$Sed_level_stand_mg <- j
    predict(modSSb4, newdata = PR_ratioSS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
PR_ratioSS_plot2 <- ggplot() +
    geom_point(data = PR_ratioSS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_reduced_PR_ratio,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Predicted probability of reduced\nP/R ratio to sediment exposure",
         color = "Binary Data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(10,max(PR_ratioSS2$Sed_level_stand_mg)), 
                  breaks=c(0.01,0.1,1,10,100,1000,loaelSS2),
                  label=c("0.01","0.1","1","10","100","1000",round(loaelSS2,digits = 1))) +
    geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelSS2, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

PR_ratioSS_plot2

By sediment exposure duration

jvalues <- with(PR_ratioSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    PR_ratioSS2$Sed_exposure_days <- j
    predict(modSSb4, newdata = PR_ratioSS2, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_dur)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
PR_ratioSS_plot3 <- ggplot() +
    geom_point(data = PR_ratioSS2, mapping = aes(
      x = Sed_exposure_days, 
      y = Binary_reduced_PR_ratio,
      color = Ref)) +
    labs(x = "Sediment exposure duration (days)", 
         y = "Predicted probability of reduced\nP/R ratio to sediment exposure",
         color = "Binary Data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_continuous(limits=c(0,max(PR_ratioSS2$Sed_exposure_days))) +
    geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability, 
                                  linetype = "solid")) +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

PR_ratioSS_plot3

Binary effects of sediment on the photosynthesis of adult corals

Lillian J. Tuttle

9/24/2020

Here I will explore how DEPOSITED AND SUSPENDED sediment affects the PHOTOSYNTHESIS of tropical, scleractinean coral ADULTS. The database to be used is comprised of data extracted from studies deemed relevant during the systematic literature review process.

The Dataset

Deposited Sediment

Suspended Sediment

Exploratory plots

DEFINITIONS:

Threshold Analyses with Binary Data

DEFINITIONS:

DEPOSITED SEDIMENT

Photosynthetic efficiency

NOAEL/LOAEL defined by exposure concentration

NOAEL/LOAEL defined by exposure duration

P/R ratio

NOAEL/LOAEL defined by exposure concentration

NOAEL/LOAEL defined by exposure duration

SUSPENDED SEDIMENT

Photosynthetic efficiency

NOAEL/LOAEL defined by exposure concentration

NOAEL/LOAEL defined by exposure duration

P/R Ratio

NOAEL/LOAEL defined by exposure concentration

NOAEL/LOAEL defined by exposure duration

Logistic Regression Model Fitting

DEPOSITED SEDIMENT

Photosynthetic efficiency

Effect estimates

Prediction figure

P/R Ratio

Effect estimates

Prediction figure

By exposure concentration

By exposure duration

SUSPENDED SEDIMENT

Photosynthetic Efficiency

Effect estimates

Prediction figures

By sediment exposure concentration

By sediment exposure duration

P/R Ratio

Effect estimates

Prediction figures

By sediment exposure concentration

By sediment exposure duration